Mergers of non-spinning black-hole binaries: Gravitational radiation characteristics 
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We present a detailed descriptive analysis of the gravitational radiation from black-hole binary 
mergers of non-spinning black holes, based on numerical simulations of systems varying from equal- 
mass to a 6:1 mass ratio. Our primary goal is to present relatively complete information about 
the waveforms, including all the leading multipolar components, to interested researchers. In our 
analysis, we pursue the simplest physical description of the dominant features in the radiation, pro- 
viding an interpretation of the waveforms in terms of an implicit rotating source. This interpretation 
applies uniformly to the full wavetrain, from inspiral through ringdown. We emphasize strong rela- 
tionships among the £ — m modes that persist through the full wavetrain. Exploring the structure 
of the waveforms in more detail, we conduct detailed analytic fitting of the late-time frequency evo- 
lution, identifying a key quantitative feature shared by the £ = m modes among all mass-ratios. We 
identify relationships, with a simple interpretation in terms of the implicit rotating source, among 
the evolution of frequency and amplitude, which hold for the late-time radiation. These detailed 
relationships provide sufficient information about the late-time radiation to yield a predictive model 
for the late-time waveforms, an alternative to the common practice of modeling by a sum of quasi- 
normal mode overtones. We demonstrate an application of this in a new effective-one-body-based 
analytic waveform model. 

PACS numbers: 04.25. Dm, 04.30.Db, 04.70.Bw, 04.80.Nn 95.30.Sf, 95.55.Ym 97.60. Lf 


I. INTRODUCTION 

The final merger of two black holes having comparable 
masses will produce an intense burst of gravitational radi- 
ation, and is expected to be one of the strongest sources 
in the gravitational wave sky. Mergers of stellar black 
holes are key targets for ground-based detectors such 
as LIGO, VIRGO, and GEO600, and knowledge of the 
merger waveforms is an important component of improv- 
ing the detectability of such systems. The space-based 
LISA detector will observe mergers of massive black holes 
at high signal-to-noise ratios, allowing a tests of General 
relativity in the strong-field, dynamical regime. 

Today, numerical relativity studies are beginning to 
progress toward a full description of black-hole binary 
merger systems. For non-eccentric inspirals, this space is 
spanned by seven parameters: the symmetric mass-ratio 
77 = m\rri 2 / (rrii+ 1712 ) 2 , and the six combined components 
of the black holes’ spin vectors. Considerable study has 
been focused on the fiducial center point of this parame- 
ter space, the case of equal-mass non-spinning black hole 
mergers. After the series of breakthroughs that ushered 
in an era of rapid progress in the field [1-4], several in- 
vestigations assessing the accuracy of the available equal 
mass waveforms and applying them to data analysis were 
conducted [5-9]. 

In this paper, we undertake a descriptive study of the 
waveforms generated in the late inspiral and merger of 
black-hole binaries for the subspace of non-spinning black 
holes, parameterized only by 77 . Our study is based on a 
series of numerical simulations, discussed in Sec. Ill, cov- 
ering at least the last > 4 orbits of nonspinning black- 


hole binary mergers with mass ratios extending to 6:1 
(77 « 0.122). Several of the simulations presented here 
have already been applied in a recent paper, focusing on 
the development of a faithful analytic waveform model 
[10]. Here we provide details of these and additional 
simulations, together with considerable analysis, focused 
on providing a qualitative and quantitative picture of 
how the waveforms from non-spinning black-hole merg- 
ers depend on 77 . Nonspinning black-hole binary merger 
waveforms were previously examined in Ref. [11], but 
our analysis is novel and complementary to that work. 
Our descriptive presentation puts emphasis on the re- 
lationships between waveforms from the different mass- 
ratio cases and different harmonic modes, with references 
to Ref. [11] where related observations have been made. 
Our approach to describing the inspiral-merger-ringdown 
transition is particularly distinct, founded in a uniform 
approach that describes all stages of this process in sim- 
ilar terms, and ultimately suggesting a complementary 
physical picture. 

Black-hole-binary merger waveforms have been noted 
for their “simplicity.” For the nonspinning systems the 
simple physics of the coalescence is exposed by a spheri- 
cal harmonic decomposition of the waveforms. In Sec. IV 
we walk readers through the basic features of the radi- 
ation, characterizing amplitude and phase evolution of 
the multipolar components, and discussing relationships 
among the simulations representing different mass ratios, 
and among the multipolar components of each simula- 
tion. As we analyze the waveforms we develop a concep- 
tual interpretation of the basic waveform features. In this 
interpretation we consider the structure of an implicit ro- 
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tating source , which could have generated the measured 
radiation through its rotational motion. This allows a 
uniform interpretation that applies throughout the coa- 
lescence process: inspiral, merger and ringdown. 

In Sec. V, we examine the strong final burst of radi- 
ation beginning ~ 20M before the formation of a com- 
mon horizon. We quantitatively describe the phasing in 
terms of an analytic model, based on a continuous, mono- 
tonically increasing frequency. We find, in particular, 
that the peak rate of change in frequency, appropriately 
scaled, is the same across all £ = to modes and mass 
ratios. We also identify relationships among the mode 
amplitudes and phases, which are connected to an ap- 
proximately linear relationship between angular momen- 
tum and frequency: cPJ/du ; 2 * * * ~ 0. We interpret these 
relationships in terms of the implicit source. 

Finally, In Sec. VI, we demonstrate the utility of what 
we have learned in our waveform characterization by ap- 
plying some of the quantitative features we have uncov- 
ered in a new variation on the analytic waveform model in 
[10], which was based on the Effective-One-Body (EOB) 
resummation of the Post-Newtonian(PN) approximation 
to inspiral dynamics [12]. In particular, we provide a dis- 
tinct late-time waveform model, alternative to the com- 
mon “spectroscopic” model [13, 14] based on sums of 
quasinormal mode overtones. 

II. OVERVIEW 

We begin with some examples of gravitational strain 
waveforms as they might be observed by gravitational 
wave instruments. In observational work, and PN analy- 
sis, it is customary to describe the radiation in terms of 
gravitational wave strain, h. In representing the strain, 
it is convenient to combine the two real waveform polar- 
ization components, h+ and h x , into one complex strain 
waveform, 

h = h+ + ih x ■ (1) 

We decompose the strain waveforms measured on a 
sphere of radius R ex t, into spin- weighted spherical har- 
monic components, he m . The details of the decomposi- 
tion, and how the waveform information is extracted from 
the numerical simulations, are given in Appendix A. 

The waveforms in this section are aligned in time and 
phase so that the maximum strain amplitude occurs at 
t = 0. The remaining figures of this paper will be aligned 
in a similar way, but with t = 0 marking the time of peak 
(2,2) mode energy flux, E 2 2 (unless stated otherwise). 

Fig. 1 shows waveforms from mergers of nonspinning 
black holes for various mass ratios, as observed at dis- 
tance R on the rotational/orbital axis of the system. 
The figure shows h + for each of the four mass ratios 1:1, 
2:1, 4:1 and 6:1. For these observers the observed wave- 
forms will be circularly polarized, so that h x is 90° out 
of phase with h + . We use units in which G = 1 and c = 1 
and express both time and spatial distances in terms of 


the total mass M, where M ~ 5 x 10 6 (M/Mq) sec ~ 
1.5(M/M 0 )km. 

More typically, the observer will not be located on the 
system’s orbital axis. The left panel of Fig. 2 shows h + 
for the 4:1 case. The strain is measured at an azimuthal 
angle of 0° and various inclinations 1 . The detailed shapes 
of the waveforms change as the system is reoriented so 
that the observer moves off the system’s rotational axis. 
For larger inclinations (closer to being viewed edge-on) 
there are notable modulations at half the base gravita- 
tional wave frequency. 

The right panel of Fig. 2 shows h + for different mass 
ratios oriented at an inclination of 90° and an azimuthal 
angle of 0°. For this orientation, h+ constitutes the full 
strain waveform. For larger mass ratios, the lower fre- 
quency modulation increases. For gravitational wave ob- 
servations of sufficiently strong binary black hole sources, 
the types of differences shown in Fig. 2 could be exploited 
to estimate the inclination and mass-ratio of the source 
system. 

For observational purposes, the combined waveform in- 
formation encoding all possible source orientations can be 
conveniently represented in terms of spin-weighted spher- 
ical harmonic components (see Eqs.(A2)-(A4)), providing 
a neat description of the leading waveform features, The 
multipolar decomposition is even more valuable as a tool 
for exposing the hallmark simplicity of the merger ra- 
diation. The readily apparent simplicity in the waves 
viewed from the system’s orbital axis in Fig. 1 extends 
to each of the spherical harmonic components. Viewed 
off axis, these components linearly combine to yield the 
more complex appearance of the waveforms in Fig. 2. 

This characterization of the gravitational radiation 
from a merging black hole binary in terms of circular 
polarization was first recognized in the Lazarus project 
studies [15]. In this picture, the radiation can be repre- 
sented by a slowly varying amplitude and a polarization 
angle; see Eq. (2) below. This description relies on how 
the radiation appears to distant observers located on the 
rotational axis of the system. Other observers will typ- 
ically see elliptically polarized waves, having a generally 
simple pattern that conforms to the rotational nature of 
the source. In the equatorial plane, the radiation reduces 
to the plus polarization, corresponding to the observer 
seeing no circulation in the source. Looking along the 
negative rotation axis, the observer sees circular polar- 
ization with the opposite helicity. 

Each of the spherical harmonic waveform components 
exhibits circular polarization with steadily varying phase 


1 The inclination angle is defined here as the angle between the 

line-of-sight with respect to the detector and the orbital axis of 

the binary. This is the same angle referred to as “inclination” in 

the PN/NR literature, and most equations are constructed using 

that definition. However, the astronomical literature has often 

defined inclination to be the angle between the line of sight and 
the orbital plane of the binary, resulting in a 90° inconsistency. 
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FIG. 1: The plus component of the strain, h+, is shown for different mass ratios using the t = 2 and £ = 3 modes. The observer 
is located at distance R along the axis of the system, at inclination 1 9 = 0° and azimuthal angle 0 = 0°. The strains are scaled 
by symmetric mass ratio and aligned such that, for each mass ratio, the peak of /122 occurs at t = 0. The phases are rotated 
such that the phases are 0° at t — 0. 



FIG. 2: Waveform variation with inclination. The left panel shows how the wave shape changes with inclination 6 for a 4:1 
mass ratio. h+ is plotted, using the sum of multipolar modes up to l = 5. The curves are labeled by the inclination in degrees 
of the observer, and the time axis is labeled so that the peak of /122 occurs at t = 0. The right panel shows how the wave 
shape changes with mass ratio. The waveforms h+ are computed at an inclination of 90° from the 1 = 2 and l = 3 modes. The 
different mass ratios are aligned in time so that the peak of /122 occurs at t = 0. The rotational phases are matched to 0° at 
t = 0. In both panels, the azimuthal angle is 0°. 


and amplitude, providing a natural framework for devel- 
oping a practical and intuitive understanding of binary 
black hole merger radiation. Our basic waveform de- 
scription, in Sec. IV, and the more detailed analysis that 
follows, is based on this spherical harmonic decomposi- 
tion. 

As we describe the waveforms, we will also suggest 
a heuristic interpretation of what the radiation tells us 


about the motion and structure of the binary black hole 
source. In the weak- field description of radiation from 
a rotating object, the multipolar waveform components 
of the gravitational radiation can be associated with dy- 
namics of multipolar moments of the radiating source 
[16]. It is useful, in conceptualizing the full coalescence 
radiation from inspiral through merger and ringdown, to 
think of the multipolar radiation description as provid- 
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ing information about the motion of a changing source 
object, described as a sum of several multipolar mass 
moments. This source object is what we will interpret as 
an effective rigid rotator radiation source, with a slowly 
changing structure. We refer to this as the implicit ro- 
tating source (IRS). 

In the process of coalescence, the source begins as a 
separated black-hole binary system and ends as a single 
distorted black hole. For nonspinning binary mergers, 
numerical and PN results consistently indicate that the 
radiation is circularly polarized, in the sense first rec- 
ognized in the Lazarus project studies [15], not only in 
the inspiral, but uniformly through the merger and ring- 
down. In our conceptual source description, this pattern 
of circular polarization is consistent with radiation gen- 
erated by rotational motion of each source multipolar 
moment, where the polarization phase is tied to the in- 
stantaneous orientation of the source. Similarly, we think 
of the amplitude of the radiation multipole as related, 
through some generalization of the quadrupole formula, 
to the magnitude of the implicit source multipole. 

We can write each multipolar component in a specific 
polar form natural for circularly polarized radiation: 


rhim(t ) 


(to > 0) , 

H lm e im <rlW (_i)< (m < 0) . 


Each amplitude H( m is expected to be a slowly varying 
function of time and can be conceptually considered as 
a function only of the magnitude of the source multipo- 
lar moment and its rotational frequency. The additional 
sign for some to < 0 cases allows a consistent interpreta- 
tion for the component phases and the component am- 
plitudes. The waveform phase is given here by m4>^, 
but in most of our analysis we refer to which we 

call the rotational phase. In our terms of our implicit 
source heuristic, the rotational phase for some particular 
(£, to) mode can be thought of as the azimuthal orienta- 
tion of the specified multipolar source component. In the 
inspiral, where the source can be considered as a sepa- 
rated binary, should coincide with the orbital phase, 
independent of (£, m) . In the post-Newtonian expansion 
[17], all the defined here agree to at least 2PN order, 
while the amplitudes Hg m remain real and non-negative. 
The equatorial-plane symmetry of these mergers ensures 
that Hi m = H^_ m ) and that <f>^ = so that we 

need consider only to > 0 modes for this analysis. 

The expansion (2) is not appropriate for the to = 0 
modes. This points to an important caveat to our implicit 
rotating source interpretation of the radiation, that it ap- 
plies to the degree that the radiation is circularly polar- 
ized. While not strictly vanishing, the m = 0 waveform 
components, and other deviations from circular polariza- 
tion are generally extremely small, and largely unmea- 
surable at the resolutions of the simulations we study. 
For the most part, we will not address deviations from 
circular polarization, and the consequent limitations of 


our implicit rotating source interpretation in this paper, 
focusing for now on the dominant features of the radia- 
tion. 


III. SIMULATIONS 

Our analysis is based on four simulations, representing 
mass-ratios 1:1, 2:1, 4:1 and 6:1. Results from the 1:1 and 
4:1 simulations have appeared in previous publications 
([5] and [10], respectively). More recently, higher accu- 
racy simulations have been presented by other groups for 
the 1:1 case [6, 18]. Our older waveform is sufficiently 
accurate for our present purpose, to examine the general 
features of the non-spinning merger waveforms. 

Our numerical simulations are carried out with 
the HAHNDOL evolution code [19], which uses finite- 
differencing methods to solve a 3+1 formulation of Ein- 
stein’s equations on a Cartesian grid. For initial data 
we solve the elliptic equation given by Brandt and 
Briigmann for conformally flat data in which the black 
holes are represented by punctures [20]. This is per- 
formed numerically using the multigrid solver AMRMG 
[21], which is second-order-accurate but tuned to give 
truncation errors typically much smaller than those pro- 
duced by the evolution code. The momentum parame- 
ters are chosen according to the 2PN-accurate quasicircu- 
lar approximation given by Kidder [22], which has been 
found to result in low eccentricity. We evolve these data 
using the moving puncture method [4, 23] with a modi- 
fied version of the BSSN equations [24, 25]. Specifically, 
as suggested in [26], we replaced the conformal factor 
variable <fi with exp(— 2 <(>), which vanishes at the punc- 
tures. Further, we added the constraint-damping terms 
suggested in [27], and the dissipation terms suggested in 
[28, 29]. For the gauge we use the specific 1+log lapse 
and Gamma-freezing shift conditions recommended for 
moving punctures in [30]. 

Accurate simulations require adequate spatial resolu- 
tion near the black holes (length scales ~ M ) as well 
as in the wave zone where the gravitational waves are 
extracted (length scales ~ (10 — 100)M). To this end, 
the grid has multiple refinement levels, determined adap- 
tively near the black holes, but fixed in regions farther 
away (typically, \x\ > 30 M) where the waves are ex- 
tracted; all grid refinement is handled within the frame- 
work of the software package PARAMESH [31]. The adap- 
tive mesh refinement criterion near the black holes is 
designed to keep the scale of the square root of an in- 
variantly defined curvature component, known as the 
Coulomb scalar [32, 33], roughly constant with respect 
to the grid spacing. Interpolation in guard-cells between 
refinement regions is fifth-order-accurate, coupling with 
differencing stencils to yield at least fourth-order accu- 
racy in the bulk. 

Spatial derivatives are taken by sixth-order-accurate 
differencing stencils, with the exception of advection 
derivatives, which are handled by fifth-order-accurate 
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mesh-adapted differencing for greater stability [34] 2 * . 
Time integration is performed with a fourth-order Runge- 
Kutta algorithm. 

The initial configurations of the simulations we analyze 
are given in Table. I. In each case, the initial separation 
was chosen to be large enough to result in at least five 
orbits. The finest resolution, hf , ranged from M/32 to 
3M/224, as required to adequately resolve the black hole 
with the smaller mass in each case. The outer boundary 
was typically at |x| > 1000M, far enough away to pre- 
vent reflections from reaching the wave-extraction region 
during the simulation. 


We have measured the individual black-hole masses 
using the apparent-horizon mass nii, a quantity calcu- 
lated from the area of each hole’s horizon, which we lo- 
cate using the AHFinderDirect code [35]. From these 
horizon masses, we calculate the symmetric mass ratio 
?y = m\m- 2 ./(rn\ + TO 2 ) 2 . This gives the most precise 
specification of the actual mass ratio attained in our sim- 
ulations. In the text we will refer to the simulations by 
the mass ratio ( e.g . 4:1). 

We define the total, infinite-separation, mass M of the 
system as an analogue for the total rest mass parameter 
used in PN studies. We measure M in two ways: 

M A H = TOi+m 2 , (3) 

the sum of the individual BH horizon masses, and 

Mrest = Madm - E b , (4) 

defined as the difference between Madm, the total energy 
of the initial data, and the (negative) binding energy of 
the binary. The binding energy E b is estimated from an 
effective-one-body PN treatment [12], given the initial 
angular momentum Jq. The result shows a very close 
correspondence between Mar and M les t , with differences 
at the level ICR 4 . For the rest of this paper, we use 
M = Mar, except for the 1:1 simulation data, where 
Mar was not available for technical reasons. 

In interpreting the late-time radiation, it is valuable 
to know the mass and spin of the final Kerr black hole 
formed by the merger. We discuss the state of the final 
black hole, determined consistently by several means, in 
Appendix C. 

For the 4:1 mass ratio case, we have carried out runs at 
three different resolutions in order to assess the quality 
of the simulations. The convergence of the constraints 
and waveforms is discussed in Appendix B. 

The most important products of our simulations are 
the gravitational radiation waveforms, which we extract 


2 Sixth-order center-differenced advection is unstable, and sixth- 

order lopsided advection is too costly in terms of PARAMESH 

guardcells, which motivated our particular modification. 
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FIG. 3: The (2,2) (top panel) and (4,4) (bottom panel) strain- 
rate waveform modes from the 4:1 mass-ratio case. Three 
resolutions are shown, time- and phase-shifted to match at 
their respective peak amplitudes. Excellent agreement is seen 
to persist throughout most of the simulation. 

from the evolved simulation data as explained in Ap- 
pendix A. Strain-rate waveforms for the 4:1 case at vari- 
ous resolutions are shown in Fig. 3, where the times and 
phases have been shifted to agree at the moment of peak 
energy flux, as is generally done in our analysis below. 
We can get some measure of the error in the waveforms by 
comparing the difference between the high and medium 
resolution simulations. Fig. 4 shows the relative differ- 
ences in amplitudes, scaled by the high-resolution result. 
Ignoring the high frequency noise, the (2,2)-mode differ- 
ences (upper panel) indicate a combination of a secular 
amplitude difference and a sinusoidal effect, which results 
from the combination of the eccentricity in the orbital 
dynamics and the difference in peak time due to limited 
resolution. These combine to give differences generally at 
the 3% level, somewhat smaller at late times. The eccen- 
tricity plays less of a role in the (4,4) differences (lower 
panel), as the relative secular error is much (~ 5 times) 
larger. Sinusoidal eccentricity effects are also visible in 
the phasing error (Fig. 5). Overall, we find waveform am- 
plitude and phase errors to be consistent with between 
fourth- and fifth-order convergence (see Appendix B). 

Assuming fourth-order convergence, and using 
Richardson extrapolation, our nominal expectation for 
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TABLE I: Physical and numerical parameters of the initial data for all the runs presented. mi, p and m 2 , p are the puncture masses 
of the two pre-merger holes, ro and Po are the initial coordinate separation and (transverse) linear momentum, respectively, 
giving rise to a total initial orbital angular momentum Jo. hf is the spatial resolution of the highest refinement level for each 
run. Madm is the total energy of the initial data. The total, infinite-separation mass M of the system is measured in two 
ways - Mah, the sum of the initial (apparent) horizon masses of the two holes, and M rea t, the sum of the ADM energy and the 
binding energy from effective-one-body theory [12]. Finally, r/ is the resulting symmetric mass ratio, as determined from the 
two holes’ horizon masses. 


mass ratio hj 

un.p rri 2 ,p ro Po Jo 

Madm Mah M res t V 

1:1 M/32 

0.4872 0.4872 10.800 0.09118 0.9847 

0.9907 1.0005 0.2500 

2:1 3M/160 

0.3202 0.6504 8.865 0.09330 0.8271 

0.9889 0.9989 0.9990 0.2228 

4:1 3M/224 

M/64 
3M/160 

0.1890 0.7900 8.470 0.06957 0.5893 
0.1890 0.7900 8.470 0.06957 0.5893 
0.1890 0.7900 8.470 0.06957 0.5893 

0.9929 1.0003 1.0004 0.1601 
0.9929 1.0003 1.0004 0.1601 
0.9930 1.0003 1.0005 0.1601 

6:1 M/64 

0.1338 0.8490 8.003 0.05559 0.4449 

0.9942 1.0000 1.0001 0.1226 
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FIG. 4: The relative amplitude error (A me dium — Afi n e)/Afi ne 
for the (2,2) (top panel) and (4,4) (bottom panel) strain-rate 
waveform from the 4:1 mass-ratio case. 


these simulations leads to an error estimate for the high- 
resolution simulation applied in our analysis of ~ 1.2x 
the difference shown in Fig. 4. To be conservative, we 
could instead assume second-order convergence, which 
would lead to an error estimate of ~ 2.8 x the difference 
shown in Fig. 4. 

The errors for the 2:1 case should be comparable to the 
4:1 case. The resolution for the 6:1, scaled by the smaller 
black hole’s mass is about 15% lower than lowest resolu- 
tion 4:1 simulation, suggesting errors eight times larger, 
if we conservatively assume fourth order convergence and 
that the errors around the smaller black hole dominate. 
The errors for the 1:1 mass ratio case are discussed in 
Ref. [5], 
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FIG. 5: The absolute phase error (^medium — <Pfine) (in radians) 
for the (2,2) (top panel) and (4,4) (bottom panel) strain-rate 
waveform from the 4:1 mass-ratio case. 


IV. DESCRIPTIVE ANALYSIS OF 
WAVEFORMS 

In this section we provide a descriptive analysis of the 
waveforms from our simulations. We try to serve two 
purposes in analyzing the radiation. In the first place, 
we are hoping to provide material for gravitational wave 
observers, and others outside the field of numerical rela- 
tivity, which makes clear some of the general character- 
istics of the radiation from these mergers. Beyond that, 
we also push the analysis in more detail, hoping to gen- 
erate deeper insight into the physics which generates the 
radiation. Through this analysis we explore the similar- 
ities and differences for the various mass-ratio simula- 
tions, and among the different multipole components of 
the 4:1 case case. In this way, we examine the waveform 
amplitudes and energy, and the waveform phasing. As 
we proceed, we will interpret the results in terms of our 
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implicit rotating source model, building up a heuristic 
description that applies through the inspiral, merger and 
ringdown of the binary. 

Following a brief discussion of strain-rate in IV A, we 
study the waveform amplitudes and the associated en- 
ergetics of the merger in IV B. In IV C we address the 
polarization phase of several waveform modes, relating 
them to a common implicit source phase. Next, in Sec. V, 
we will examine the late-time frequency evolution and the 
relation to amplitude in more quantitative detail through 
the inspiral-merger-ringdown transition. 


A. Strain-rate 
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In the introduction, we motivated the spherical har- 
monic phase and amplitude waveform decomposition 
with a discussion in terms of strain h. In analyzing our 
numerical simulation results, however, we can work more 
directly with the strain-rate h(t) = dh(t)/dt (see Ap- 
pendix A for a more detailed discussion). As with the 
strain decomposition (2), we will expand the strain-rate 
as 


Rh/tm (t) 


(to > 0) 

-iA im e irn ' s '^ (- 1 )* (to < 0 ) , 


with A (rn real and non-negative. 3 Direct differentia- 
tion of (2) reveals the relationships between the phases 
and amplitudes defined in (2) and (5), with A^ m ~ 
+ 0(5PN) while differs from <I+ m only 
at 2.5PN order. Note that the differentiation produces a 
phase shift of 7t/2, so that strain-rate phases should be 
defined in the slightly unconventional form (5) if we wish 
to preserve the property that all are equal to orbital 
phase in the limit of well-separated binaries. This rep- 
resentation allows a more meaningful comparison of the 
phases of different multipolar modes, ensuring that in- 
stantaneous phase corresponds to the orientation of the 
binary system during the inspiral. 

The strain-rate amplitude is directly related to the ra- 
diation power for the (£, to) mode by Eq. (A7), E# m = 
(A^ m ) 2 /167r. Henceforth we shall use the strain-rate- 
derived rotational phase &e m , rather than We also 
limit our presentation to the to > 0 modes, as equatorial 
symmetry implies = $q_ m ) and A lm = A^_ m y 

Unless otherwise indicated, in the remainder of the pa- 
per, the time axis of each plot will be shifted so that the 
peak of E 22 (and hence of the strain-rate amplitude A 22 ) 
occurs at t = 0. As the (2,2) mode is strongly dominant, 
this will closely approximate the peak time of the total 
E. 


FIG. 6: The (2,2) mode of the peak power versus mass-ratio 
T] (solid circles), with a numerical fit to the quadratic-quartic 
formula (6). Also shown is the full peak power, divided by 
two (asterisks). The growing difference between the two as 
77 —> 0 indicates the increasing importance of the subdominant 
modes in unequal-mass cases. The right-hand scale gives the 
peak power in c.g.s units. 


B. Amplitude and energetics 

We first study wave amplitudes across modes and mass 
ratios. Since we are examining strain rate waveforms, the 
modal energy flux is effectively equivalent to the square of 
the mode amplitude, as in Eq. (A7). Preferring the most 
physical language we will express our modal amplitude 
comparisons as energy flux comparisons. In terms of the 
implicit rotating source model, we can think of the energy 
carried by the radiation as energy lost by the source. 

In Fig. 6, we plot the actual peak values of the dom- 
inant (2,2) energy flux contribution from Eq. (A7) as 
a function of symmetric mass ratio r). As this mode 
contribution is proportional to I/122I 2 , and we expect 
lim^o 1 1222 1 — » 0 (the test-particle limit), we fit it to 
a quadratic-quartic form, obtaining the fit 4 : 

E 22 W) = (4.40 ± 0.17) x KTV 

+(5.43 ± 0.31) x 10“ V. (6) 

We also plot the peak of the total energy flux for each 
mass ratio, scaled by one-half, since the (2,2) and (2,-2) 
modes contribute equally to E. The difference between 
Atotal/2 and E 2 2 increases as 77 —> 0, reflecting the 
increased importance of other modes for unequal masses. 

Aside from the value of the radiation maxima, it is in- 
teresting to see how the radiation power evolves in time 
near the peak. In Fig. 7, we show shapes of the domi- 
nant (2,2) contributions to the peak energy fluxes (A7) 


3 A similar expression, differing from (2) by an overall sign, would 

be equally appropriate for direct interpretation of numerically 4 A quadratic-cubic form is equally plausible, but fit the numerical 
derived ij > 4 waveforms. data worse in this case. 




for each mass ratio. These are scaled to the same peak 
height to allow shape comparison, and shifted in time so 
that the peaks of E 22 are aligned. We note the strik- 
ing similarity of the peak shape and duration across all 
mass ratios. During the late-inspiral phase, the more ex- 
treme mass ratios appear to radiate more energy; how- 
ever, since we have normalized each curve by peak height, 
this only means that the equal-mass binary experiences 
a steeper climb to its peak power rate. Nevertheless, the 
different mass ratios follow similarly shaped tracks ap- 
proaching merger, and differences have been effaced by 
~ 10M before the peak power. The post-peak portion 
of the curve is determined by the dominant quasi-normal 
mode (QNM) damping time, which varies only slightly 
with underlying Kerr spin for moderate spins (see, for 
instance, tables in [13]). The lower spins on the black 
holes formed by smaller-?; mergers should cause power to 
fall off faster in these cases. The inset in Fig. 7 shows 
the difference in fall-off rate relative to the equal-mass 
case from 20 M after the peak. 

In Fig. 8, we concentrate on the 4:1 mass ratio, plot- 
ting the strongest modal contributions to the flux, again 
scaled to the same peak height. Here we have applied 
only an overall time shift, so that the total flux peak is 
at t = 0. The energy profiles of the radiation burst in 
modes with t = m are similar, all peaking at approxi- 
mately the same time. The subdominant modes are rel- 
atively stronger in the burst than in the inspiral, so that 
they show up here as weaker in the approach to the peak 
even after scaling to the peak radiation. We note that 
the subdominant modes with l = m, (3, 3) and (4, 4), are 
particularly similar in this regard. Similarity in £ = m 
modes, and distinction in the other modes is a general 
feature of the bursts in several ways. 

The shape of the peaks with / / m are particularly 
distinct. The (2,1) mode peaks particularly late, and the 
burst is much stronger than the inspiral. We note that 
the (3,2) mode shows a double-bump in its contribution 
to the energy flux. From Fig. 9, this appears to be robust 
in its gross shape over resolution and extraction radius. 
From Fig. 10 we note, however, that the extent of this 
double-bump effect is very dependent on mass ratio; it 
is less in the 6:1 case, and not evident at all in the 2:1 
case. Later we will also note irregularities in the late-time 
frequency evolution of this mode, apparently indicating 
a deviation from circular polarization in this case. 

In both Fig. 7 and Fig. 8 we have normalized the 
mode-flux peaks for the purposes of shape comparison. 
It is also important to understand the relative strengths 
of each mode. We show in Fig. 11 the relative mode con- 
tributions to E(t) for several dominant modes over the 
final inspiral and merger of the 4:1 case. In the merger- 
ringdown peaks, as in the inspirals, the £ = m modes 
dominate the energy flux, followed by the £ = |m| + 1 
modes. More discussion of the relative mode strengths 
for general nonspinning mergers is given in Ref. [11]. 

In addition to the direct mode contributions, we plot 
the PN-derived energy “partitions” - the power emitted 



FIG. 7: The shape of the (2,2) contribution to E(t) for all 
mass ratios, scaled to unity at the peaks. We have shifted the 
time axis so that the peak in 15(2,2) occurs at t = 0. The inset 
focuses on the ringdown portion of the curves. 



FIG. 8: The shape of the different multipolar contributions 
to E(t) for the 4:1 case, scaled to unity at the peaks. We 
have shifted the time axis so that the peak in 15(2,2) occurs at 
ip 0. 

in each mode as a fraction of the total, according to the 
leading-order “restricted” PN expressions found in Eqs. 
(30)-(36) of [10], where the underlying orbital frequency 
was derived from the (2,2) mode. These are shown in Fig. 
11 using dashed lines. Until near the peak time, where we 
discontinue the PN curves, we see that the partitioning 
tracks the numerics very well for all modes except (2,1), 
which grows visibly faster than the restricted amplitude 
prediction after ~ 100M before merger. A similar study 
[10] showed that the restricted (leading-order) approxi- 
mation for the amplitudes consistently overestimates the 
strength of the radiation. This shows that the leading 
order partitioning of energy can provide a simple, but 
more accurate, approximation of the mode amplitudes. 
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FIG. 9: The (3,2)-mode “double-bump” in E(t) for the 4:1 
case. We show the inner two extraction radii (20A7 and 45 M) 
for each of our three resolutions, shifting the time axis so that 
the peak in E( 2 , 2 ) occurs at t = 0. Though there is variation 
in the detailed shape, the overall double-bump envelope seems 
to be robust. 



FIG. 10: The (3,2)-mode in E(t) for the 2:1 (red/dashed), 4:1 
(green/dot-dashed), and 6:1 (blue/dotted) cases, extracted at 
Rext = 45 M. We have shifted the time axis of each data set so 
that the the peak in U( 2 , 2 ) occurs at t = 0. The amplitudes are 
unsealed. It is clear that the 4:1 case has the most pronounced 
deviation from a single well-defined peak. 

We will take advantage of this in Sec. VI to provide more 
accurate amplitudes in a variation on the analytic EOB- 
based waveform model studied in [10]. 


C. Waveform phasing 

In gravitational-wave observations, the waveform 
phase provides most of the time variation in the signals, 
and consequently is critically important in encoding ob- 


FIG. 11: Relative power Ee :7n for different multipolar modes 
of the 4:1 case, on a logarithmic scale (solid lines). For com- 
parison, we also plot (dashed lines) the fractional power ex- 
pected from each of these modes based on Eq. (A6), and the 
restricted PN waveforms of Eqs. (30)-(36) of [10]. The differ- 
ent modes have been aligned to peak at t = 0. The right-hand 
scale gives the power in e.g.s units. 

servable information about the source. Here we consider 
the phasing of the leading spherical harmonic waveform 
components. Following the discussion above, we concep- 
tually interpret each waveform phase as describing the 
orientation of a particular (£, m) multipole of an implicit 
rotating source of the gravitational waves. 

Direct comparative analysis of phases provides a 
stronger probe of the phase relations among the multi- 
polar modes than the comparative analysis of frequencies 
conducted in a number of previous studies of numerical 
simulations. Here we will discuss waveform phasing in 
terms of < f>^ m as it appears in (5), which we will compute 
from each strain-rate mode h( m (t). As noted in Section 
II, we expect all phases to agree in the large-separation 
limit. 

We first compute the strain-rate waveform phase ipe m 
using the conventional decomposition: 

Rh l m{t) = V im {t)e i ^ t \ (7) 

with V( m real and non-negative. Then, setting this equal 
to (5), and solving for the rotational phase we 

find 

$ = f m (Wm(t) - f + 2tt n tm ) (to > 0) 

tm 1 m ivirn{t) + § + 27T7 l£ m + tl t) (m < 0) . 

( 8 ) 

The ![ term results from the factor i in (5), while the 
(— l) 1 factor there produces the lit term for to < 0. The 
ne m terms express the 2-7T ambiguity in defining wave- 
form phase. Considering any {£, to) mode in isolation 
leads to an m-fold degeneracy in the associated rota- 
tional phase. We resolve this degeneracy by choosing the 
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pair {n22,ri33} £ {(0,1) x (0,1,2)} that yields the clos- 
est consistency between $22 and $33 at early times (near 
t = —400 M). We then determine the remaining n^ m for 
closest early consistency with 4>22- This gives us a phase 
for each mode that can be interpreted as the rotational 
phase of the implicit rotating source that produced that 
component of the radiation. 

Fig. 12 shows the rotational phase from several 
modes of the strain-rate waveforms from the highest- 
resolution (3M/224) 4:1 run, together with the rotational 
phase calculated from the tracks of the punctures. The 
left panel shows that for the inspiral portion of the evo- 
lution, all waveform phases, extracted at i? ext = 45 M, 
agree extremely well, except for the (2, 1) and (3, 2) 
modes, which differ by a significant part of a radian. The 
relative difference of each mode from the (2, 2) mode is 
shown in the right panel. The differences between the 
t = m modes are < 7 t/ 60, much smaller than the 27r/m 
ambiguity in defining rotational phase from the wave- 
forms. 

Note that the early part of the (2, 1) waveforms con- 
tain, by far, the longest wavelength radiation present, 
suggesting a greater potential for problems caused by 
extracting the waveforms too close to the source, not 
yet in the wave zone. For the relatively deviant (2, 1) 
and (3, 2) modes, we performed Richardson extrapolation 
with respect to extraction radius R ex t, using the values 
extracted at 45M and 90 M and assuming an R~^ t error 
for each mode. These Richardson-extrapolated phases 
are also shown in the right panel of Fig. 12, subtracted 
from the (2,2) rotational phase. Richardson extrapola- 
tion evidently reduces the early phase deviations in these 
modes considerably. 

The rotational phases calculated from the £ = m modes 
agree to within 0.025 radians during the inspiral for sev- 
eral hundred M before the peak E, effectively identical 
within the uncertainties of the numerical approach. This 
phase agreement is consistent with expectations based on 
the PN analysis, for which all (£, m) phases should agree 
up to 2PN order. These rotational phases also agree to 
within about 0.2 radians with the coordinate-dependent 
rotational phases measured from the puncture tracks, af- 
ter shifting the puncture track phase by an overall factor 
of 7t/2. Heuristically, we can think of each (£,m) ra- 
diation waveform mode as having been generated by the 
rotation of its own implicit source component. The phase 
agreement would then be interpreted as indicating that 
these implicit source components remain aligned through 
the inspiral. This is to be expected for a system which 
can be effectively described as an orbiting pair of point 
particles. 

It is, perhaps, more remarkable that a very tight agree- 
ment among the £ = m mode persists throughout the 
merger and even into the ringdown, remaining within 
about 1 radian until ~ 50M after the merger, when the 
amplitude has already diminished significantly. Accord- 
ing to our interpretation, this phase agreement suggests 
that a significant portion of the implicit rotating radia- 


tion source maintains some structural integrity through- 
out the coalescence. That is, the implicit source we 
have considered appears to exhibit considerable “rigid- 
ity” through merger. This is only possible because of the 
close relationship among the fundamental £ = m quasi- 
normal ringdown frequencies (see Sec. VA), mimicking 
the harmonic frequency relationship that holds during 
the inspiral. For the £ ^ m modes this quasinormal 
frequency relation does not hold and the phases must 
separate in the merger. In terms of our implicit source 
picture, these £ m components of the source seem to 
shear away from the main source structure to rotate at a 
faster rate. 

Heuristically, the puncture motion is strongly tied to 
the rotation of the implicit source for most of the evo- 
lution. During this period, it is natural to think of the 
implicit source as an inspiralling pair of point-like ob- 
jects moving on time-like world lines. At late times the 
orientation phase angle of the puncture track disassoci- 
ates from the waveform rotational phase. The punctures 
veer away from the implicit source at a late times as they 
fall into the final black hole. At this point, though we 
can continue to consider an implicit rotating radiation 
source, it no longer makes sense to think of that source 
as a pair of point-like objects. 

Having compared the phases of different multipolar 
modes for the 4:1 case, we now consider how the phase 
evolution depends on mass-ratio. There are various rea- 
sonable approaches to comparing the phases among sim- 
ulations of the different mass-ratio cases. Having estab- 
lished above the rotational phase consistency for the dif- 
ferent (£, m) modes, we will compare only the dominant 
(2, 2) phases. An obvious approach is to compare phases 
directly against time, scaled by the total PN mass M. In 
the early-time well-separated limit, however, the leading- 
order PN analysis indicates that phases for different mass 
ratios should evolve at similar rates when time is scaled 
by the chirp mass M. = M77 3 / 5 

Fig. 13 shows the rotational phase computed from the 
(2,2) mode of strain-rate for different mass ratios. In the 
left panel, we align the rotational phases at an early time 
and scale time by the chirp mass M. For this plot, we 
shift the rotational phases in (chirp) time so that at t = 0, 
the chirp frequency, which is rotational frequency multi- 
plied by chirp mass, is 0.033 and the rotational phase 
is 0. Following this approach, we would expect good 
phase agreement at times sufficiently early that only the 
leading-order PN effects are significant. However, for 
the late portion of coalescence that we have simulated, 
we find that the different mass ratios remain roughly in 
phase for several hundred M before and after t = 0, 
peeling away in order at late times. 

In the right panel we compare phases in a manner com- 
mon for numerical relativity waveform comparisons, we 
shift the curves in time so that each peak energy flux 
occurs at t = 0, and we rotate the phases so that the 
phases are 0 at this time. We scale the time by M. For 
the equal-mass case Ai = 0.435275, and for the other 
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FIG. 12: Rotational phase calculated from puncture tracks and different multipolar strain-rate modes extracted at 45 M 
for a 4:1 mass ratio. The rotational phase is computed from each mode using (8) and from the angle of the vector connecting 
the two punctures. (The puncture tracks remain in the x-y plane.) The phases are aligned so that the peak of E occurs at 
t = 0 for all phases. The left panel compares the different calculations of while the right panel compares the differences 
between the phases and the rotational phase from the (2,2) mode . In the right panel, we add Richardson extrapolations of the 
(2,1) and (3,2) modes based on the 45 M and 90 M extractions of these modes and assuming an R~ 2 error. The Richardson 
extrapolations are distinguished by asterisks in the figure. The differences are smoothed to better show the trends. 



FIG. 13: Rotational phase $22 for different mass ratios, as computed from /i 22 - In the left panel, we time shifted the curves 
such that the chirp frequencies Q 22 M = Q 22 M are equal at t = 0. In the right panel, we aligned the curves such that the /122 
peaks occur at t = 0, and the phases are 0 at this time. 


mass ratios it is smaller, so all of the curves in the left 
panel are stretched by at least a factor of 1/M ~ 2.3 in 
time relative to the curves in the right panel. In the M- 
scaled right panel, the different mass ratios again remain 
approximately in phase for several hundred M before and 
after t = 0. At sufficiently late times, and particularly 
for small 77 , we might expect this manner of consistent 
phasing as the evolution of the system eventually (after 
the innermost stable circular orbit (ISCO)) may become 
dominated by the course of unstable geodesic trajectories 
around the larger black hole (or an effective black hole in 
the effective-one-body (EOB) framework). In that case 
the frequency evolves independently of the more strongly 


independent rate of energy or angular momentum loss. 


V. DETAILED LATE TIME ANALYSIS 

In Sec. IV we have presented general information about 
the phasing and amplitudes of the radiation components. 
Our analysis has stayed close to the standard numerical 
relativity waveform analysis, though we have emphasized 
an interpretation in terms of an implicit rotating source 
model. In this section we go beyond the standard wave- 
form presentation, exploring the radiation with the hope 
of developing a deeper understanding of the simple char- 
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FIG. 14: Rotational frequency for several (£, m) modes of the 
4:1 mass ratio simulation. These are time derivatives of the 
curves in Fig. 12. For £ = m the rotational frequencies remain 
similar from inspiral through ringdown. These are similar to 
the puncture track frequency until t ~ —40 M, while the (2, 1) 
mode begins to deviate somewhat earlier. 

acteristics of the radiation as described above. Those fea- 
tures and our heuristic interpretation suggest a new ap- 
proach to examining the structure of the late-time phas- 
ing and apparent relationships between frequency and 
amplitude evolution. 

In Sec. V A, we examine the phasing again, seeking 
a quantitative understanding of the the late-time evolu- 
tion of the polarization frequency. We introduce a prac- 
tical model that captures the merger-ringdown transition 
without the need for multiple quasi-normal mode over- 
tones. We investigate the implications of this model in 
relating the frequency and amplitude close to merger in 
Sec. VB. 

This section is more technical than Sec. IV, with some 
subtle discussion of late-time radiation characteristics. 
For readers who may wish to jump ahead to Sec. VI, 
we note two results that we will carry forward: (1) a 
simple quantification of the peak chirp rate to for t = m 
modes, and (2) the idea that dJ/dto becomes approxi- 
mately constant at late times, which may serve as sum- 
mary of relationships between frequency and amplitude 
near the radiation peak. 


A. Waveform frequency evolution 

First we will study the phasing in the merger and ring- 
down in more quantitative detail by comparing the po- 
larization frequency evolution for each mode with a sim- 
ple empirical model. Based on our heuristic model, we 
interpret the polarization frequency as corresponding to 
the rotational frequency (£, to) component of the implicit 
source. 

Fig. 14 shows the evolution of the rotational frequency 


for several (£, rri) modes of the 4:1 mass-ratio case. Sim- 
ilar frequency evolution curves are a common feature in 
papers on numerical relativity waveforms (see, for exam- 
ple, Refs. [1, 7, 11]). These curves are time derivatives of 
the phase evolution curves shown in Fig. 12, zoomed in 
on the late-time behavior, near the elbows in the phase 
curves. The striking similarity in phasing for the various 
t = to modes implies similar frequency evolution, which 
has been noted in previous studies [36] . 

At late times, this similarity in frequency is made 
possible because of a special approximate relationship 
among the fundamental £ = m quasinormal ringing fre- 
quencies, that they are nearly equal after dividing by 
the azimuthal mode number to to get what we call the 
rotational frequency. This has been considered in [37], 
which pointed to a connection between the quasinormal 
ringing frequencies and the frequencies of stable null or- 
bits of a black hole at the “light-ring” with frequency 
12 lr = l/(o + M(r + /M)( 3 / 2 ^). The association extends 
to charged Kerr-Newman black holes and has been com- 
pared with recent precise quasinormal ringing frequency 
calculations in [13] and [38]. Conceptually, this allows us 
to think of the rotational frequency of the l = m modes 
at late times as corresponding to the rotational rate of 
gravitational perturbations orbiting at the “light-ring” . 
This suggests a heuristic description of our implicit ro- 
tating source at late times as a gravitational distortion 
of the forming final black hole which predominantly re- 
volves around the black hole on null orbits at the light 
ring. 

Returning to Fig. 14, we note that the (3, 2) mode is 
different from all the others, showing two spikes, near 
t = 9 M and 37 M. Comparisons of waveforms extracted 
at different radii, and from simulations of different reso- 
lutions, suggest some sensitivity to extraction radius, but 
do not suggest that the features will vanish in more ac- 
curate simulations or with more distant wave extraction. 
These anomalies may be related to the unusual shape 
in the amplitude peaks noted in Fig. 9 above. We will 
discuss this mode’s behavior further in Sec. VII. 

For all other modes the frequency evolution follows a 
simple smoothly evolving curve, qualitatively similar in 
each (£, m) and mass-ratio case (see Fig. 15 below). In 
particular, we note that, except for small noise contri- 
butions, each curve shows that the frequency increases 
monotonically, ultimately saturating at a frequency set 
by the fundamental quasinormal ringdown mode. This 
monotonic frequency development is a universal charac- 
teristic of the radiation from inspiral, through merger, 
and up to ringdown. In the PN analysis of quasi-circular 
inspiral, this characteristic makes it possible to describe 
the changing structure of the hardening binary as a func- 
tion of frequency instead of the more coordinate-specific 
separation. This allows us, for instance, to write the 
waveform amplitude as a function of frequency. 

Subsequently we will assume monotonic frequency de- 
velopment throughout the coalescence process. This 
principle underlies our empirical curve fitting of the fre- 
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quency evolution, allowing more quantitative analysis of 
the late-time phasing evolution. In Sec. VB we will fur- 
ther apply this idea as we study relationships between 
late-time frequency and amplitude evolution. 

To produce an empirical curve for describing the late- 
time frequency evolution, we assume that each {£, m) 
mode has a monotonically increasing polarization fre- 
quency, which approaches the fundamental ringdown fre- 
quency wqnm at late times. The general expectation that 
the frequency decays exponentially toward the ringdown 
frequency suggests that we model frequency evolution 
based on the hyperbolic tangent function. 

Specifically, we will compare frequency evolution of the 
strain-rate waveforms with a model of the form f l(t) = 
g(t), where 


g(t) — fii + (f2f 


a) 


1 + tanh[ln + (t — t 0 ) /b] 


(9) 

This provides a curve that first grows exponentially, with 
e- folding time 6/ (2k), from some initial frequency fl,, 
then decays exponentially, with e-folding time 6/2, to 
the final frequency fif. The presence of the exponent 
re allows the early exponential growth rate to differ from 
that at late times. The early part of this model is a coarse 
approximation to the growth in frequency near the end 
of the inspiral. This approximation must, therefore, fail 
to fit the data if we look back sufficiently long before the 
time of peak radiation, but, as we will show, it provides a 
fair approximation of the approach to the peak. The rate 
at which the frequency grows, the chirp-rate, increases to 
some maximum, then decreases to zero on approach to 
the final frequency Of. For a more meaningful parame- 
terization, we set 


7 / -j \ 1 + Av 

f]; = f]f — -fi° f 1 + — j . (10) 

With this choice, dg/dt will peak with value flo at time 
t 0 . The model then depends on the five parameters re, b, 
O 0 , to and Of. 

As shown in Figs. 15 and 16, we find that fits to the 
model g(t) provide an excellent approximation to the nu- 
merical data for the strain-rate rotational frequencies in 
most significant cases. For the (2, 2) modes of all mass- 
ratios, and for all but one of the significant (£, to) modes 
for the representative 4:1 mass-ratio case, we find agree- 
ment within a few percent after t = —20 M, with the 
primary differences coming from apparent noise in the 
numerical simulations, within the uncertainties in the nu- 
merical results. 

Fig. 15 shows the comparison for the (2,2) mode ro- 
tational frequencies for several mass ratios. A glance at 
the curves shows that the unequal-mass cases are quite 
similar to the equal-mass frequency evolution, previously 
examined in Refs. [7, 39]. The dominant difference for 
the unequal-mass cases is that the final frequency fif de- 
creases with 77, consistent with the decrease in the spin 


of the final black hole produced. We expect flf to cor- 
respond to wqnm/to, where wqnm is the fundamental 
(n = 0) quasinormal ringing frequency for the specific 
(£, m) mode for a black hole with the appropriate spin. 
In the infinite-mass-ratio limit (rj — ► 0), should corre- 
spond to half the Schwarzschild (2, 2) quasinormal mode 
frequency Mwqnm = 0.3737, indicated by the horizontal 
dashed line in the left panel of Fig. 15. 

We show a few examples of polarization frequency 
curves, for subdominant modes in Fig. 16. The (3, 3) 
mode is very similar to the (2, 2) modes shown in Fig. 15, 
as are the other £ = to modes (not shown). The (2,1) 
mode is of similar shape, also well approximated by our 
fit. 

The quantitative fit results are summarized in Table 
II. The error bars are based on statistical fit estimates, 
also incorporating the ranges of best fit results obtained 
by varying the fit range starting between t/M = —25 and 
— 15 and ending between t/M = 20 and 60. The final fre- 
quencies flf approached in the fit curves in Figs. 15 and 
16 were robustly determined by the fits within a fraction 
of a percent. The fif frequencies from the fits in Fig. 15 
were applied in Table V to find final black-hole parame- 
ters consistent with those determined by conservation of 
energy and angular momentum. 

The peak in the chirp-rate flo is a particularly sig- 
nificant quantity in determining the shape of curves of 
our general form g{t). As shown in Table II, our fits 
determine flo to within a few percent in all cases. An 
interesting relationship is apparent among the fits for all 
the l = to cases for all values of g studied. In the last 
column of the table, we show the peak chirp-rate values 
scaled by the mass of the final black hole Mf and the 
final frequency flf. In each case where £ = to we find 
CloMf/Slf ~ 0.02, consistent within the fit uncertainties. 
This scaling makes some sense, since the height of the rise 
in frequency through the final radiation burst is largely 
determined by flf, while the time scale over which this 
rise occurs seems to be similar when time is scaled by 
Mf. In Sec. VI we will use this result to predict the 
phase evolution in an analytic waveform model. 

Our model for late-time frequency evolution (9) de- 
scribes exponential decay toward Of at an e-folding rate 
given by half our fitting parameter b. For all cases, the 
values of b are within about 30% of 10M. In some cases, 
the fits for b are rather sensitive to the initial starting 
time, varying by up to 20% or 30% in the (2,1) and 
(4,4) modes of the 4:1 case. At this coarse level, we note 
that the values for b are similar to the exponential decay 
rates for quasinormal ringdown mode amplitudes listed 
in Table IV. We will consider this relationship further in 
Sec. VB below. 

The other parameters in our fit are re, relating to the 
shape of our fit curve at early times, and to, giving the 
time at which the frequency peak occurs. The parameter 
re is not very precisely determined; as we would expect, it 
depends sensitively on the starting time of the fit inter- 
val, since the early exponential frequency growth is only 


14 



0.015 


0.01 


0.005 


a 

<1 


-0.005 



FIG. 15: Strain-rate rotational frequency evolution in the merger for (2,2) modes of several mass ratios (dashed lines), with an 
analytic fit (solid lines). The fit encodes a monotonically increasing frequency, which at late times decays exponentially toward 
the fundamental ringdown mode. The horizontal dashed line marks M flf = 0.18685, half the (2,2) mode QNM frequency for a 
nonspinning perturbed hole. The right panel shows the residuals, which are comparable to uncertainties in the numerical data 
after 20A-f before the peak. 



FIG. 16: Strain-rate frequency evolution in the 4:1 mass-ratio merger for several spherical harmonic modes, with an analytic 
fit as in Fig. 15. The left panel shows the polarization frequency uu m - The right panel shows the residuals of the unsealed 
rotational frequencies = utm/m. 


TABLE II: Results of fitting Eq. (9). All quantities are scaled by the final mass Mi = Mf ;I . a< j defined in Eq. (Cl). The 
parameter flf can be related to the fundamental real QNM frequency, uqnm, and has been used to extract two different 
estimates of the final dimensionless spin a of the hole - see Table V. 


mass ratio 

(e,m) 

K, 

b/Mf 

TYi£)jqM{^ 

mllfMf 

to/M 


1:1 

(2,2) 

.7 ± .1 

13.2 ± .3 

0.0112 ± .0001 

0.528 ± .001 

-4.6 ± .2 

0.0210 ± .0002 

2:1 

(2,2) 

0.6 ± .1 

12. 3 ±.8 

0.0104 ± .0002 

0.5023 ± .0006 

-3.6 ±.2 

0.0207 ±.0004 

4:1 

(2,2) 

0.5 ± .1 

10.5 ± .2 

0.0096 ± .0001 

0.4566 ± .0001 

-4.8 ± .3 

0.0210 ± .0002 


(2,1) 

1.0 ±.5 

12 ±2 

0.0120 ±.0005 

0.421 ± .006 

-3.9 ±.8 

0.0285 ±.0012 


(3,3) 

0.3 ± .1 

9.7 ± .2 

0.0153 ±.0001 

0.730 ± .001 

-3.3 ±.l 

0.0210 ±.0002 


(4,4) 

0.15 ± .1 

7.5 ± .7 

0.0212 ± .0003 

0.991 ± .005 

-1 ± 1 

0.0213 ± .003 

6:1 

(2,2) 

0.5 ± .1 

10.4 ± .5 

0.0089 ± .0002 

0.4349 ± .0004 

-6.5 ±.3 

0.0205 ± .004 
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a coarse approximation of the expected behavior. The 
values for to show that the peaks in 17 generally occur 
roughly 4 M before the total energy peaks at t pea k- As 
was the case for the power peaks in Fig. 8, the chirp-rates 
of the different spherical harmonic modes peak at slightly 
different times. 

We have supplemented our general implicit rotating 
source picture with the additional idea that the rota- 
tional frequency for each mode grows monotonically, not 
only in the inspiral, but also through the merger and ring- 
down. Based on this expectation we have identified an 
analytic fit model for the late-time frequency evolution 
that precisely matches the data for all cases but the (3,2) 
mode. These fits provide a quantitative understanding of 
the late-time phasing yielding, in particular, a robust re- 
sult for the peak chirp rate fio for all £ = to modes. We 
will apply this information in Sec. VI. 


B. Late-time frequency and amplitude 
relationships. 

The last step in our waveform analysis is to consider 
relationships between the frequency evolution and the 
waveform amplitude. 

In the PN description of the quasicircular inspiral, the 
orbital frequency not only tells us the rotational rate, 
but can also serve as a label for describing the momen- 
tary state of the rotating object (in the inspiral case this 
means that we can reference the state of the system in 
terms of r(I7)). The PN generalization of the quadrupole 
formula, describing radiation from the rotating system, 
then leads to an expression for amplitude as a function of 
frequency. Our description of the gravitational radiation 
suggests an implicit source rotating with monotonically 
increasing frequency as it continues to “harden”, as the 
system evolves smoothly into merger and ringdown. In 
this section we seek to further unify this picture of the full 
coalescence process, considering an analogue of the PN 
description of amplitude as a function of frequency that 
can describe the radiation in the merger and ringdown. 

In Sec. IV B we emphasized that the radiation power 
E( m provides essentially the same information as the 
strain-rate amplitude (A6). If the wave frequency is 
known, then the modal contribution to the total radia- 
tive angular momentum can similarly provide informa- 
tion about the gravitational wave amplitude. Eq. (A9) 
gives an expression for angular momentum flux in terms 
of wave amplitude and phase. The relation simplifies to 

Jem ~ 167r ^ (Ai m ) 2 (11) 

up to 5PN order, indicating that the radiation carries 
maximal angular momentum J = E/i 7, as we gener- 
ally expect for circularly polarized radiation. If we know 
how the rotational frequency evolves, we can derive the 
mode amplitudes Ai m from the mode-by-mode relation- 


ship of either energy or angular momentum with fre- 
quency, Ee m (i 7) or J im (t 7). 

To approach an understanding of the late-time rela- 
tionships between amplitude and frequency, we compare 
how the system’s energy and angular momentum ap- 
proach their final values, with how the system’s frequency 
approaches its final value. In Fig. 17 we examine the 
radiative loss of energy dE and angular momentum dJ 
as the coalescing system approaches its final quiescent 
state, comparing these with the remaining difference of 
gravitational- wave frequency from its late-time limit S7f , 
as determined in Table II. As well as our standard final 
rotational frequency f7 = 4> defined in (5), we also show 
the frequency: based on the strain, 17/j = 4>/i as defined 
in (2). We have also rescaled the energy and angular mo- 
mentum by a constant, selected so that the value matches 
that of f7/j at the time of peak radiation power. 

Fig. 17 indicates a general correspondence between 
how the angular momentum approaches its final state, 
and how the gravitational wave frequency approaches its 
final state. Of the several cases of spherical harmonic 
modes and mass ratios that we have examined in this 
manner, we show two examples: the (2,2) case for the 
equal-mass simulation, where the evolution of angular 
momentum and frequency correspond most closely (left 
panel), and the (2, 1) mode from the 4:1 mass ratio sim- 
ulation, with the weakest correspondence (aside from the 
nonconforming (3,2) case) (right panel). The correspon- 
dence is closest, holding to a better approximation over 
a longer period of time, in associating dJ with 17^. The 
association with energy is slightly weaker. 

Though these plots suffer significantly from small mod- 
ulations in frequency that we have not resolved numer- 
ically, the results suggest an approximate relationship 
between rotational frequency and angular momentum 
in particular. If the frequency evolution is otherwise 
known, then the late-time evolution of angular momen- 
tum for each mode could be approximately described by 
J(t) = Jf + £(I7(f) — I7f), where £ is a case-dependent con- 
stant that we will not attempt to specify generically. As 
£ = dJ/dfl we may refer to it as the dynamical moment 
of inertia of the implicit rotating source. This late-time 
expectation creates the possibility of extending our PN- 
based understanding of angular momentum flux into the 
late time waveforms, giving us the additional information 
we need for a full (approximate) description of J(f7). 

This relation between frequency, angular momentum, 
and amplitude provides a connection between the peaks 
in modal radiation power, shown in Fig. 11, and the peaks 
in chirp rate given in Table II. Assuming dJ/d^l = £, a 
constant, yields 

E « £1717 (12) 

=> E « £(m7 + I2 2 ) (13) 

Since the 17(7) curve is steeply increasing, we would ex- 
pect the peak in E to be near the peak in 17 but slightly 
delayed. Expanding (13) linearly about ten when 17 = 0, 
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FIG. 17: Relationship between angular momentum and frequency evolution approaching ringdown for the (2,2) mode of the 1:1 
case (left panel) and the (2,1) mode of the 4:1 case (right panel). We show the differences dE and dj from the final mass and 
angular momentum as measured from the radiation, together with the difference from the final ringdown rotational frequency, 
measured from h (dfi) and from the strain h (dfi^)- We see similar evolution for the dominant (2,2) mode of each mass ratio 
simulation, and for all modes other than (3,2). In the left panel the angular momentum and energy (in geometrical units) have 
been scaled up by factors of 2.5 and 10, respectively, to roughly match the frequency residuals at t = 0. In the right panel the 
rescaling factors are 760 for angular momentum and 2000 for energy. The agreement of the curves, with these rescalings, is 
consistent with dJii/ d£l 22 ~ 0.40 from t = 20 through the peak radiation. The more approximate agreement in the right panel 
is consistent with dJ 2 i/dQ- 2 i = 0.0013, through the peak radiation. 


we find that E reaches its peak at a time 

Cl 2 

tpeak ~ to — > to, (14) 

“olio 

where Ho and Ho are the frequency and its third time 
derivative at to- These can be evaluated from (9) and 
the fit parameters k, b, H 0 , and fif in Table II; we find 
that the values for t pea k — to all lie in the approximate 
range (2.16,3.17), in rough agreement with the direct fit 
for to in Table II. 

Conceptually, the monotonic evolution of frequency 
might lead us to hypothesize a relationship between the 
structure of the implicit rotating source and its rotational 
frequency, such that the changes in the structure of this 
radiating rotator will be associated with finite changes 
in frequency. Knowing that the frequency growth must 
be limited by the quasinormal ringing frequency implies a 
peak in the chirp rate. This leads to the expectation that 
finite changes in energy and angular momentum are as- 
sociated with finite changes in frequency, so that dE /dfi 
and dJ/dfl approach constants once the evolution in fre- 
quency slows. The peak in radiation power might then be 
viewed as a consequence of the peak in chirp-rate. More 
quantitatively we find that dJ/dQ seems to be roughly 
constant even before the evolution in frequency slows 
down (i.e. before t = 0). This correspondence will be 
applied in the next section to provide a model for ampli- 
tude evolution through the peak, based on information 
about the frequency evolution. 

If we could postulate the constancy of dJ/dH, we might 
also apply that assumption to “explain” some of what we 


have seen above. In Sec. V A we noted a rough agreement 
between the timescale b in our frequency fitting curve (9) 
and the quasinormal ringdown amplitude decay rates for 
the corresponding quasinormal modes. Following the dis- 
cussion above, this relationship could be derived, in the 
H — > wqnm/*w limit, from the constancy of either energy 
or angular momentum losses with respect to change in 
frequency. For constant £ = dJ/dH, Eqs. (A6), (5) and 
(12) imply that 

A 2 = \Rh\ 2 = 16? tE ss 16<HH. (15) 

In the H — > fif limit, the strain-rate amplitude decays 
at the rate predicted from black hole perturbation theory, 
A tm -> A 0 im exp (-t/Te m ), where Ti m is the e-folding 
rate for the amplitude decay for the fundamental (£, m) 
quasinormal mode. In this limit, our frequency evolution 
fit model reduces to 

H(t) fif - (H f - fi;)e- 2( *-‘ o) / 6 . (16) 

Applying these limiting expressions for amplitude and 
frequency in (15) in the limit H — > fif yields 

(A 0 ) 2 e~ 2t/T w 32 ^ f (Hf - fii)e- 2(t -‘ o) / b , (17) 

b 

where the left-hand-side derives from the amplitude, and 
the right-hand-side from frequency. Ignoring the con- 
stant coefficients, this implies that b = r. 

If we adventurously assume the constancy of £ on ap- 
proaching the ringdown, and expand the amplitude in 
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powers of e = exp (— t/b ) = exp (— f/r), the implied am- 
plitude frequency relation might also provide more infor- 
mation about the amplitude evolution. Since D(t), and 
consequently the right-hand-side of (17), contains only 
even powers of e, the next term in the expansion for am- 
plitude should be 0(e 3 ). This suggestion motivated our 
expansion (C4) applied in fitting the late-time amplitudes 
in Sec. III. 


VI. VARIATIONS ON THE EOB MODEL 

An approach to modeling black-hole binary radiation 
known as the effective-one-body (EOB) model has been 
presented in the literature [12, 40-45]. The late-time 
waveforms in these models are based on a now-common 
description of the merger process as an epoch of radia- 
tion from spiraling particle-like trajectories, followed, in 
a sudden transition, by black hole ringdown dynamics 
with waveforms described by a superposition of quasi- 
normal frequencies. Our waveform analysis provides a 
complementary description of black-hole binary merger 
radiation that can be applied in an alternative late-time 
waveform model. 

A recent promising approach along these lines is the so- 
called pseudo~4PN (p4PN) EOB model [10]. This model 
extends the 3PN-accurate EOB metric with a term of 
4PN order with a tunable multiplier A. The phasing 
obtained from this expansion is combined with leading- 
order PN strain amplitudes (the restricted approxima- 
tion) to obtain mode-by-mode waveforms valid for inspi- 
ral. For the merger and ringdown, the phasing and ampli- 
tude are derived simultaneously from a superposition of 
quasinormal modes. For each {£, to) angular mode, the 
fundamental ringdown modes and a few overtones are 
summed in proportions as required for continuity with 
the late end of the radiation from the inspiral phase. The 
value of the p4PN multiplier A is then chosen to match 
the pre- and post-merger waveform portions, optimizing 
the agreement with full- numeric waveforms. 

In this section, we show that it is possible to develop 
variations on the p4PN EOB model that usefully encode 
some of the waveform phase and amplitude relationships 
we have described above. A key difference with the new 
variant is our prescription for the transition from inspiral 
to merger-ringdown radiation. In contrast with [10], for 
each (£, to) mode we consider the entire wavetrain as that 
of a slowly varying instantaneously rigid rotator, consis- 
tent with the dominant “circular-polarization” waveform 
pattern encoded in the radiation. The phase evolution 
might be thought of as arising from the rotation rate of 
the corresponding (£,m) source structure, which is con- 
tinued through ringdown by continuous matching to a 
function of the form (9). The wave amplitude will be de- 
rived directly from expectations for the energy or angular 
momentum content of the radiation. 

We present two specific models of this nature. Model 1 
is based on exactly the same EOB-based prescription for 


inspiral-plunge trajectories in [10], while Model 2 shows 
the effect of a slight variation in the underlying EOB 
model. In both variants, as in [10], we derive the wave- 
form phasing directly from the EOB trajectories (with 
A = 60 for the strength of the p4PN term) up until some 
matching time, which we take simply as the time at which 
the (2,2) wave frequency is half the frequency of the fun- 
damental (2,2) ringdown mode. After this point we use 
our fit model (9) to describe the subsequent phase evo- 
lution. 

Recall that this model depends on several parameters: 
k, b, Do, to, and Df. The results of our analysis in Sec- 
tion V B guide us in producing a fully specified model for 
these parameters. We take Df = wqnm/w from the fun- 
damental ringdown frequency wqnm of the radiation. For 
the time-constant for frequency decay b we use the fun- 
damental quasinormal mode amplitude decay time con- 
stant. While this is not clearly implied by our fits in the 
last section, it will lead to the correct amplitude fall-off, 
as specified below. For the strongest l = m modes, our 
fits indicate Do = 0.021ujQNM/wAff . Lacking any bet- 
ter model we simply increase this by a factor 4/3 when 
£ ^ to, roughly consistent with the higher value of Do 
found for the (2,1) mode. For these models we derive 
the quasinormal modes using the fit for the final black 
hole mass and spin described in [10]. The remaining pa- 
rameters to and k are chosen to provide continuity, up 
to the second time-derivative of phase, with the direct 
EOB-based phasing at the matching time. 

Our prescription for the wave amplitudes differs from 
the restricted amplitude description applied in [10]. As 
we showed in Fig. 11, we can improve on the restricted 
amplitude approximation by using the leading order PN 
expressions for waveform (JL, to) mode amplitudes only 
to fix the partitioning of radiation power into angular 
modes. We set the total power independently, from the 
full-order EOB model description of the radiation power; 
the resulting waveforms are then energetically consistent 
with the EOB description of the dynamics, and also show 
better agreement with the numerical results (note that 
in this model there is no radiation in the nonrotational 
m = 0 modes). After reaching the matching frequency, 
we continue the amplitude evolution based on the as- 
sumption, suggested in Section VB, that the amplitude 
is roughly constant through the radiation peak. 

In our new model, we set the late-time amplitude by 
asserting the approximate relationship (15), written this 
time in terms of the polarization frequency ui = toD: 

a 2 ^ -i a r /i o\ 

A(,m ~ 16 TT&m ^2 j (18) 

setting the value of for amplitude continuity at the 
match-frequency. With this model for the amplitude, 
the peak in the gravitational wave amplitude is a direct 
consequence of the peak in the time-derivative of gravi- 
tational frequency, fixed by Do in (9). The exponential 
decay in amplitude also follows directly from the expo- 
nential approach of the wave frequency to oaqnm- 
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FIG. 18: Comparisons of new EOB-based models with numerical results for the (2,2) waveform frequency and amplitude for 
the 1:1 (left) and 4:1 (right) mass ratio cases. The models are variations on the p4PN EOB waveform model with a flux-based 
determination of the wave amplitudes and an alternative, based on Eq. 9, to quasinormal mode summing for continuing the 
waveforms through ringdown. The vertical bar indicates our matching frequency where we transition from direct EOB phasing 
and flux. Models 1 and 2 correspond to different versions of the radiative flux, which primarily affects the consequent wave 
amplitudes; see the text for more details. 




FIG. 19: Comparisons of new EOB-based models with numerical results for the 4:1 case, for the (3,3) (left) and (2,1) (right) 
waveform modes. The plots show frequency and amplitude, labeled as in Fig. 18 


We compare the frequency and amplitude of the mod- 
eled waveform (2,2) component with the corresponding 
numerical result for the 1:1 and 4:1 mass ratio cases in 
Fig. 18. A similar comparison is shown for some of the 
subdominant modes in the 4:1 case, in Fig. 19. The model 
we have so far described, based on the p4PN EOB tra- 
jectories is labeled Model 1 in the figures. The matching 
frequency is indicated as a vertical line in each plot. The 
frequency curves indicate very good phasing agreement 
for all cases except the (2, 1) mode of the 4:1 mass ratio 
case. In that case the sharp rise in frequency occurs a 
few M too late, several times worse than the agreement 
shown for the other modes. The difference is a conse- 


quence of the slightly higher frequency early on, as the 
numerical (2, 1) frequency begins to grow already before 
the matching point as compared with the rotational fre- 
quency consistent with the numerical £ = m modes. 

The generally good amplitude agreement shown be- 
fore merger represents an improvement over the simple 
restricted amplitude model employed in [10]. For the 
£ = to modes, Model 1 overestimates the amplitudes just 
before and consequently after matching. This indicates 
that the PN-based flux expression applied to the p4PN 
EOB model, based on [46] and [47], overestimates the 
flux at high frequencies. This is perhaps not surprising, 
since the flux formula is formulated with a pole at finite 
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frequency. Physically we expect the flux to decrease at 
late times, when the frequency approaches the quasinor- 
mal ringdown frequency. 

To correct for this, we show also results for a second 
variation, Model 2, in which we have introduced a zero 
in the flux function at ringdown frequency. Following the 
notation of [46] , we modify their Eq. (45) for the flux to 
read 


Fpn 


32 

5 


rj 2 v w 


1 — w/wqnm 

1 - v/l’pole 


fp N {v\ri), 


(19) 


and then respecify the coefficients C\ — C7 in their Eq. 
(50) to again provide consistency to 3.5PN order with 
the Taylor series expansion for the flux. Note that the 
flux now depends on £ and to via the quasinormal mode 
frequency wqnm and v = (Mfl) 1 / 3 = (Mw/m) 1//3 . The 
modified flux function anticipates that the radiation will 
cut off at the ringdown frequency. We find that the new 
EOB model provides a very good approximation to the 
original p4PN EOB phasing (with A = 60) if we choose 
A = 27 for the new version with modified flux. The 
figures show that, for the £ = m waveform modes, the 
amplitudes based on Model 2 with the modified flux show 
better agreement with the numerical results leading up 
to the matching point, and also at the peak, while the 
frequency evolution is nearly identical. For the (2,1) case, 
Model 2 suffers the same problems as Model 1. 

We have applied the observations made in previous 
sections to successfully predict late-time waveforms. This 
provides an alternative description of the transition to 
ringdown, distinct from the widely-applied approach of 
summing quasinormal modes [7, 11, 48]. In Ref. [10], we 
applied it to nonspinning mergers. 

This raises the question: how is it that each of these 
quite distinct approaches to approximately describing the 
late-time radiation can be simultaneously effective? We 
can explore this question by considering the polarization 
frequency evolution in a waveform constructed from a 
sum of quasinormal modes (suppressing £ and to labels) 


s{t) = J2 A n e~ ant+lV,n(t) , (20) 

n— 0 

with ip n (t) — ip n ( 0) + ui n t and where uj n and a n = l/r n 
correspond to the n th quasinormal overtone mode. Next 
we restrict to the first two terms and to linear order in 
e(f) = exp[— (<Ji — <To)f], which vanishes at late times. 
This yields 

s{t) « A 0 e~ a « 4 +^o(t) ^ 1 + ee i(vi(t)-vo(t))^ (21) 

~ Aq e -CT ot+e c °s(Avj)+i[<po+esin(Avj)] (22) 

where Act = ui — tio and A ip(t) = — <po(t). Taking 

the derivative of the expression in square brackets gives 
the polarization frequency of s(f) 

u} s (t ) « coo — e [Act sin (A ip) — Aw cos (A<^)] , (23) 


where Aw = wi— wo- Note that the expression in brackets 
is periodic with period 27 t/Aw, where Aw is the differ- 
ence between the fundamental quasinormal ringing fre- 
quency and its first overtone (see [13] for a table of these 
overtones) . 

For the waveform modes we consider here MAw is 
quite small, generally 1 or 2%, which means that the pe- 
riod of the expression in brackets is > 100M. We rewrite 
w s , only keeping terms linear in Aw, as 


w s (f) 


Wq — eAcrsin (A(^) 



Aw 

Act 


cot {Aip) 


(24) 


Considering the second term, we note that MAa is gen- 
erally just under 0.2 for the cases we’ve studied so that 
(Aw/Acr) < 0.1. The other factor, cot (Aip), is not pre- 
dictable without knowledge of the initial conditions, but 
most values of Aip leave the second term somewhat less 
than one. 

If these conditions hold near the onset of ringdown, we 
might neglect the time dependence in the second term. 
The late-time frequency evolution is approximated by a 
exponential decay with a time constant of l/Aer. Work- 
ing under these assumptions, and comparing with the 
late-time frequency evolution model in this paper, yields 
the association 2b = l/Aer, where b is the fitting param- 
eter in our late time frequency evolution model. Again, 
looking at the quasinormal mode values, we comment 
that eq/cr 0 ~ 3 so that l/(2Acr) ~ b ~ r 0 , the approx- 
imate relationship noted in Sec.V A, which was also ap- 
proximately derived from the assumption that dJ/dio is 
constant at late times. Consistently, applying the rela- 
tion t7i/cro ~ 3 in Eq. (22) also give an expression for 
late-time amplitude consistent with Eq. (C4). 


VII. DISCUSSION 

With hope of reaching out to a wide range of re- 
searchers interested in gravitational radiation from black- 
hole binary mergers, we have provided a descriptive walk- 
through of many of the general features of late-time wave- 
forms from generic mergers of nonspinning binary black 
hole systems, based on a series of numerical simulations 
covering systems from equal-mass up to mass ratio 6:1. In 
this basic waveform description we have examined wave- 
form phase and amplitude, comparing results among dif- 
ferent mass ratios, as well as among the different spin- 
weighted spherical harmonic (£. to) component modes. 

In our presentation, we have attempted to describe 
the radiation in the simplest physical terms, pointing out 
traits in the waveforms that are similar through the in- 
spiral, merger and ringdown stages. Throughout the co- 
alescence, we find simple waveforms in each (£, to) mode, 
each exhibiting strong circular polarization and mono- 
tonically increasing polarization frequency. 

In our amplitude comparisons, we find that the leading 
order PN-prediction for energy-partitioning provides a 
good estimate of the amplitude until late in the merger 
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for £ = m modes. In astrophysical units, our fit for the 
peak-power in the dominant (2, 2) mode is 

E 2 Av) = (1-6077 2 + 19.70?7 4 ) x 10 67 ergs -1 . (25) 

Scaling the amplitudes by the peak values yields a very 
similar shape through the peak for the (2, 2) mode ampli- 
tudes for all mass-ratios we have studied. For a particular 
mass-ratio case, the peak-widths remain similar, though 
there is some variation in peak-time. 

For each mass ratio, the phase (and frequency) of the 
different (£, to) components are strongly related. While 
this should be expected for the early inspiral, where the 
waveform phase is directly connected to the orbital phase, 
we show that, for l = m modes, the same relationship 
holds through the merger and into the inspiral. We com- 
pare the phasing among simulations with different mass 
ratios in two ways, with time scaled by chirp-mass, as is 
appropriate in the early-time limit, and then with time 
scaled by total system mass M. With the latter scaling, 
the dominant (2, 2) waveforms are similar in phase (and 
77-scaled amplitude) for the last ~ 200M. 

In the near-peak waveform comparisons the (3, 2) 
mode does not exhibit the same simple behavior as the 
other modes. This mode is easily subject to coupling 
with the much stronger (2, 2) mode. It has been seen 
in [7, 36] that the (3,2) mode demonstrates significant 
mode mixing with the dominant (2,2) mode - the QNM 
ringdown part of the (3,2) waveform contains the funda- 
mental frequencies of both the (3,2) and the (2,2) modes. 
We speculate that this mixing might be partly due to the 
use of coordinate extraction spheres that are systemat- 
ically warped from the areal-radius spheres appropriate 
for correct radiation extraction. A more refined choice 
of extraction spheres and perhaps a better tuned decom- 
position basis (e.g. spheroidal harmonics) would make 
it possible to represent the (3,2) waveform content in a 
manner which is like that seen with the other modes. 
Thereby we expect that a similar simple physical repre- 
sentation of the waveform content can be extended to all 
in / 0 modes. 

We suggest a simple conceptual interpretation that ap- 
plies through the full coalescence. We think of the radi- 
ation as being generated by an implicit rotating source, 
with each [£, to) mode generated separately by the {£, to) 
moment of some implicit source (which we understand 
here only in the context of the radiation). The nearly 
fixed relationship among the [£, to) phase moments is 
interpreted to indicate that the implicit source main- 
tains some structural integrity throughout the coales- 
cence, without shearing among the various modal compo- 
nents. For the £ = m modes, this rigidity is maintained 
through the merger and into the ringdown, a relation- 
ship made possible by the approximate equality for each 
of the £ = to quasinormal modes wqnm/w ~ UlRj where 
I2 lr is the orbital frequency of unstable circular prograde 
graviton (or photon) orbits. 

The following physical picture may underlie these rela- 
tionships. For well-separated black holes, the fields that 


embody the implicit source object evidenced in the ra- 
diation may be tied directly to the point-like centers of 
the orbiting black holes. The source rotation frequency is 
the orbital frequency of the time-like trajectories traced 
out by the black holes. As the binaries spiral together, 
the pair can continuously be viewed as a shrinking, dis- 
tributed dumbbell-like rotator. Eventually, most of this 
dumbbell shrinks inside the light-ring, which roughly co- 
incides with the potential barrier in the wave mechanics 
of gravitational perturbation theory. From inside, little 
radiation can escape to a distant observer, and the time- 
like motion of the black hole centers disconnects from the 
radiation. At late times, the effective radiation source be- 
comes a gravitational disturbance orbiting the forming 
black hole at the light ring. This is a seamless transi- 
tion, with nearly consistent rotational phasing among all 
£ = to modes throughout the process. For the £ 7^ m 
modes the associated quasinormal-ringing dynamics are 
somewhat distinct, and the phasing and amplitude rela- 
tionships begin to peel away from the main £ = m trend 
through the merger process. 

For the late-time portions of the waveforms, includ- 
ing the approach to the peak and the ringdown, we have 
introduced a quantitative fitting model based on a mono- 
tonically increasing polarization frequency for each mode, 
which decays exponentially toward the expected funda- 
mental quasinormal ringdown frequency wqnm at late 
times. These fits provide an excellent match for the 
frequency evolution beginning ~ 20M before the peak, 
and allow precise estimates of wqnmj as well as the peak 
rate of change in frequency ujq. Scaling the latter quan- 
tity by the final black hole mass Mf and wqnm, we find 
woAff/wQNM ~ 0.021 for all £ = to waveforms we have 
looked at, including the (2,2) modes of each mass-ratio, 
and modes up to (4,4) for the 4:1 mass-ratio case. 

Conceptually, the monotonicity of the frequency evolu- 
tion suggests that, as is the case for inspiralling systems, 
the frequency can be taken to label the state of the adi- 
abatically changing implicit rotating source that we in- 
terpret as the source of the radiation. If we suppose that 
changes in the source structure are tied to loss of angu- 
lar momentum, then we would expect that finite changes 
in frequency would be associated with finite angular mo- 
mentum loss, so that dJ(u)/dhj has a finite, non-zero 
value even at late times. Since u> approaches a non-zero 
constant at late times, we would likewise expect dJ/doj 
to approach a constant value. In Sec. VB, we show that 
the late-time evolution of J and to are approximately con- 
sistent, mode-by-mode, with constant dJ/dto beginning 
about ~ 20 M before merger. 

Such a relation between frequency and angular mo- 
mentum also implies a connection between frequency and 
amplitude. The moment of peak amplitude is expected 
to be near the peak in uj(t). At very late times, constant 
dJ I die implies a connection between the rates at which 
the frequency and amplitude approach their quasinormal 
late time state, namely that our fitting parameter b ss t, 
where r is the damping time of the quasinormal ampli- 


21 


tude decay (see Eq. 9). 

The simple relationships between the waveform modes 
and simple dependence on mass-ratio make it possible to 
specify much of the late-time waveform information de- 
veloped in our numerical simulations in terms of just a 
few quantities. This information can then be combined 
with information from the PN approximation about the 
inspiral trajectories to provide analytical models for full- 
coalescence waveforms. In Sec. VI we have applied the 
PN-consistent EOB-based p4PN trajectory model pre- 
sented in [ 10 ] together with assumptions asserting sev- 
eral of the approximate waveform features observed in 
Sec. IV. As in [10], early waveform phasing is derived 
from the p4PN EOB trajectories, with waveform am- 
plitudes based on our PN-basecl power partitioning, to- 
gether with the PN flux model. As the waveform ap- 
proaches the anticipated peak, we match to a waveform 
phasing model based on our fit model, with parame- 
ters specified according to the approximate relationships 
identified in Sec. IV, and with amplitudes derived from 
dJ/dw=constant. The fits show excellent phase agree- 
ment with the numerical simulation results for the most 
significant £ = m modes. A variation in the PN flux 
model that enforces that the flux for each (£, m) mode 
vanish as w — * ojqnm gives better late-time amplitude 
agreement. 

Lastly, we observe that our description of the late- 
time phasing and amplitudes provides a picture comple- 
mentary to another approach applied in several previous 
studies [7, 10, 11, 48], which successfully treat the late- 
time waveforms as a sum of quasinormal fundamental 
and overtone modes for each {£, m) waveform compo- 
nent. This is motivated by the expectation that wave- 
forms from generic initially compact distortions of the 
forming black hole will quickly reduce to a sum of these 
quasinormal harmonics [49]. In [10], we have shown that 
in comparisons with some of the runs presented here, 
for mass ratios up to 2 : 1 , this assumption can lead to 
a predictive waveform model with similar accuracy to 
our alternative model presented in Sec. VI. As a link 
between the two approaches, we have shown that at late 
times the combination of the fundamental and first QNM 
overtones for a particular mode may, under reasonable 
circumstances, mimic the amplitude decay properties of 
our model. 

This work suggests several directions for further study. 
In the immediate future, we plan to assess the fidelity 
of the available nonspinning waveforms and models, and 
the impact of mass ratio on the overall detectability of 
the merger signal. We also plan to apply our implicit- 
rotating-source description as a baseline in analyzing fu- 
ture higher-precision numerical simulations. This might 
provide insight into understanding finer features of the 
merger physics, some of which could violate our simpli- 
fied description. Further understanding the anomalous 
(3, 2) mode waveforms will be a first step in this direc- 
tion. We must also investigate whether our description of 
the merger radiation applies also to spinning black holes. 


It is plausible that even precessing systems might be an- 
alyzed in this way using a spherical harmonic basis that 
tracks the orbital axis [50]. This may make it possible 
to extend the analytic EOB-based waveform model pre- 
sented here to include spin effects. Including spins in 
such analytic models will be necessary for observational 
data analysis applications. 
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APPENDIX A: RADIATION EXTRACTION 

We extracted radiation using the outgoing Weyl scalar 
ip 4 , defined as in [51], calculated with a symmetric tetrad. 
This is related to the complex gravitational wave-strain 
h{t, x) via 

ip 4 (t,x) = —h + (t,x) + ih x (t,x). (Al) 

ip 4 , the strain h , and its time-derivative h (which we 
call the strain-rate) are all functions of time t, extrac- 
tion radius R ex t, and polar angles 9 , </>. As is customary 
in numerical relativity, we decompose the radiation into 
spin-(- 2 )-weighted spherical harmonic components: 


OO t 


E E C tm (t,R)- a YT(0,<f>). 

m——£ 

(A2) 

00 £ 

E E h em (t,R) -2Y e m (e,<p), 

£=2 m——£ 

(A3) 

00 £ 

E E h tm(t,R) -2Y e m (e,<P). 

1-1 m——£ 

(A4) 


With this, two time integrations of our measured quan- 
tity, ip^, yields the more familiar gravitational wave strain 

h. 

To extract the radiation from a simulation, we define a 
series of coordinate spheres of different radii i? ext ; here we 
use extraction spheres having radii between R ext = 40A7 
and R ext = 100 M. We extracted the radiation in modes 
by integrating ip 4 against different - 2 Y™( 9 , (p) over these 
coordinate spheres, using fourth-order interpolation onto 
each sphere followed by Newton-Cotes angular integra- 
tion. 
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The gravitational waves produced by the binary carry 
both energy and angular momentum. Overall, the rate 
of energy emission is given by an angular integral of the 
squared strain-rate \h\ 2 over a coordinate sphere (see Eq. 
(5.1) of [15]): 

— — = lim (f dil\h\ 2 . (A5) 

dt oo 16? r / 1 1 v ’ 

Then using Eq. (A3), we can express the total energy 
flux (A5) as a sum over modes: 

Rf/m = 


where we have used the strain-rate decomposition (5) and 
taken the limit R — > oo to go from (A6) to (A7). 

Similarly, the rate of radiation of the ^-component 
of angular momentum can be expressed as a sum over 
modes [52]: 

/ d J Z \ 777 

Jem = ( —rr ) = lim — R 2 Im(h e . rn h* f , m ) . (A8) 

V dt / l m R — > 00 lb 71 " 

Substituting the definitions (5) and (2) into expression 
(A8) for the angular momentum, and taking the limit 
R — > oo yields 

jpm = Yr COs(??r(^^ m (A9) 


= i im 

dt ) lm R-^ 00 lo7T 

(A, m ) 2 

167T 


(A6) 

(A7) 


APPENDIX B: CONVERGENCE 

We carried out three runs of the 4:1 mass ratio model 
at different resolutions to study the convergence proper- 
ties of our simulations. For these cases, the mesh spacing 
of the finest grids (the ones including the smaller punc- 
ture) was taken to be hf = 3M/160 (low resolution), 
hf = M/64: (medium resolution), and hf = 3M/224 
(high resolution). To facilitate comparisons among these 
cases, the overall grid structure of the runs was kept the 
same. In this Appendix, we discuss the convergence prop- 
erties of the constraints and gravitational waveforms. 

In comparing our medium and high resolutions for the 
4:1 mass ratio case, the Hamiltonian constraint was found 
to be manifestly fourth-order-convergent in the dynam- 
ical strong-field, where the black holes move, evidently 
dominated by the expected fourth-order error from re- 
finement interfaces. The convergence falls off to an ap- 
parent rate closer to first order in the coarsest regions. 
This seems to result from stronger dissipation of high 
frequency noise. Though less noise is generated in the 
higher resolution simulations, a greater portion of it sur- 
vives propagation into the distant coarser regions. As 
shown below, this does not appear to affect the wave- 
forms, which are well-resolved in the wave zone. The mo- 
mentum constraint appeared to be at least second-order- 
convergent in the dynamical strong-field region, but also 



FIG. 20: Convergence of strain-rate amplitude for (2,2) (up- 
per panel) and (4,4) (lower panel) modes. The medium-high 
differences are scaled for both fourth- and fifth-order conver- 
gence. Three-level convergence is lost around 150M before 
merger (vertical dashed line). 


fell off to an apparent rate closer to first order in the 
coarsest regions - from the wave-extraction region out- 
wards. 

In our waveform mode analysis, we have typically time- 
and phase-shifted the data so that the amplitude peak 
fell at time t = 0 [1]. In Fig. 3, the real part of the (2,2) 
and (4,4) strain-rate harmonics are shown for all three 
resolutions. The agreement between the resolutions is 
then seen to be excellent. In Figs. 4 and 5, we showed 
the errors to be expected from the maquillaged strain- 
rate waveform data. 

However, the presence of considerable eccentricity in 
the binary makes is difficult to compare such time-shifted 
waveforms between resolutions and establish an unam- 
biguous order of convergence. In Fig. 20 we show the 
amplitude errors for our three resolutions without time- 
shifting - that is, we plot the data starting from the ini- 
tial time in each case, scaling the medium - high differ- 
ences for fourth and fifth-order convergence. It is clear 
that we observe convergence between fourth and fifth 
order throughout the evolution until around 150M be- 
fore merger, when the difference in merger times among 
the runs becomes important. Additionally, it is possible 
that at this higher- frequency stage of merger, the lowest- 
resolution data is no longer in the convergence regime. 
We do note, however, that the rate of growth of ampli- 
tude difference between the medium and high-resolution 
runs is comparable here to that observed closer to the 
merger for time-shifted data (see Fig. 4). From this we 
deduce that the medium and high-resolution runs are still 
in the convergence regime, with errors consistent with 
fourth-order convergence. 
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FIG. 21: Convergence of strain-rate phase for (2,2) (upper 
panel) and (4,4) (lower panel) modes. The medium-high dif- 
ferences are scaled for both fourth- and fifth-order conver- 
gence. Three-level convergence is lost around 150A7 before 
merger (vertical dashed line). 


APPENDIX C: END STATES 

In this appendix we discuss the state of the final black 
hole formed in our simulations. This can be measured by 
several independent means, which we compare, finding 
agreement to within 0.4% for Mf and 2.1% for a. Table 
III, we present “coarse” results for the simulations. The 
total radiated energy Ai? ra( j is obtained by integrating 
Eq. (A5), and the total radiated angular momentum 
A J ra( j is obtained by summing over mode contributions 
(A8). The final mass and angular momentum of the post- 
merger black hole are calculated using 

Mf jrad = MadM — AUrad, (Cl) 

df,rad = JO Ao/fad, (C2) 

where, by the symmetries of the current simulations, we 
only deal with the z component of angular momentum. 
For the 1:1 simulations, we only had the leading-order 
(2, ±2) radiation modes available, so our A J rad estimate 
is significantly truncated, by as much as 11% (a con- 
servative error estimate based on the effect of similarly 
truncating mode contributions past £ = 2 for the 4:1 
case); we have marked this and derived values. We also 
quote the measured value of f pea k, the time at which E 
reaches its peak. All waveform and derived plots in this 
paper have been time-shifted by subtracting this time, 
as an approximate marker of the time of merger, unless 
otherwise indicated. 

Table III also contains data about the common ap- 
parent horizon (CAH) of the merged binary: the time 
fcAH at which the CAH was first detected, the CAH’s 
irreducible mass Mi rrj cAH, and its full (horizon) mass, 


obtained from the Christodoulou [53] formula: 

J 2 

Aff 2 CAH = M 2 rI + 2 , (C3) 

^iiK/irr 

where we use J = Jf, r ad for the final hole’s angular mo- 
mentum. 

We can use the data from Table III to estimate the 
mass and dimensionless spin of the end-state Kerr black 
hole through different methods. Given the final mass es- 
timates A/f. rad and Mf.cAH, we can calculate a = Jf/Mf. 

Another means of characterizing the final black hole 
comes from studying the characteristics of the radiation 
after the peak. This is expected to be a sum of the 
hole’s quasi-normal modes (QNMs). We can determine 
the imaginary part of the mode frequency by fitting the 
waveform amplitude Ai m from (7) post-peak to a de- 
caying exponential. The damping time in this fit corre- 
sponds to the imaginary frequency. Specifically, we fit to 
a functional form: 

A{t) = A 0 e" t/rQNM (l - C'e" 2t/TQNM ) . (C4) 

The presence of the additional damping term, 
parametrized by C, allows more freedom for non- 
linear decay early in the ringdown. Note that this factor 
of three between primary and secondary damping times 
approximately mirrors the difference in damping times 
one sees between the first two QNMs of a Kerr hole 
[54]. The results of this fit are given in Table IV. The 
value of t qnm> together with an estimate of the real 
QNM frequency wqnm, uniquely determines the mass 
and dimensionless spin of the Kerr hole. We combine 
our present fit for tqnm with the real frequency wqnm, 
as presented in Table II, to obtain (Mqnm>&qnm)- 

If we already have a robust determination of the Kerr 
mass, then we only need one component of the QNM 
frequency to lock down a. In Table II, we use wqnm to 
obtain a QNM ,re, assuming M qnm = M f)rad . 

We have gathered all the estimates of Mf and/or a 
discussed above in Table V. It is instructive to compare 
these values with recent predictions for the final mass 
and spin from analytic and numerical methods. Working 
with a subset of the data supplied here, and also invoking 
the test-mass limit, [10] suggest the one-parameter fits: 

Mf/M = 1 - (1 - y/8/9)r) - (0.498 ± 0.027)?7 2 ,(C5) 
a = \/l2r7 - (2.900 ± 0.065)r7 2 . (C6) 

The numerical data from Table V also fits well the 
wholly numerical formula (3.17a) from [11], as well as 
the parametrized formula of [55] . 

We note that when discussing the mass Mf of the end- 
state Kerr hole in Sections V A and later, we have con- 
sistently taken the radiation-derived value, Mf = Mf |rad , 
as results from [5] indicate this has errors up to a factor 
of three less than QNM-based mass. Since a far larger 
proportion of the initial angular momentum than mass is 
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TABLE III: Results of radiation and apparent-horizon analysis of merger. The total energy flux E (A5) reaches its peak at 
time fpeak/AL AE Ta d and AJ ra d are the total energy and (z) angular momentum radiated during the simulation (the latter 
calculated as a sum over modes (A8)), resulting in final mass Mf jrad and angular momentum Jf jlad (For the 1:1 case, only 
the 1 = 2 modes were available, so the emitted AJ ra d will underestimate, and Jf jlad will overestimate, the physical results.). 
Mj rrj cAH is the irreducible mass of the common apparent horizon, first detected at simulation time tcAH- A/^cah is the mass 
deduced from this and Jq rad using the Christodoulou formula [53]. 


mass 

ratio 

h f 

Rext/M 

tpeak/M AR rad (x 10 _2 )/M Aff, rad /M 

A J rad (x 10 _2 )/M 1 2 3 4 5 6 7 8 9 10 J f , rad /A/ 2 

tcAn/M M irr ,CAH/M M[ : CAh/M 

1:1 

M/32 

60 

1303.7 

3.5934 

0.9548 

33.79 f 

0.6468 f 

1196 

— 

— 

2:1 

3M/160 

45 

627.7 

2.8306 

0.9606 

23.34 

0.5937 

581 

0.9063 

0.9637 

4:1 

3M/224 

45 

641.9 

1.4327 

0.9786 

12.78 

0.4615 

588 

0.9489 

0.9796 


M/ 64 

45 

652.6 

1.4262 

0.9786 

12.99 

0.4594 

599 

0.9489 

0.9793 


3M/160 

45 


1.4102 

0.9789 

12.90 

0.4603 

610 

0.9492 

0.9797 

6:1 

M/ 64 

45 

564.9 

0.9212 

0.9850 

07.83 

0.3666 

513 

0.9667 

0.9851 


TABLE IV: Fit results for tqnm from Eq. (C4). tqnm is 
used for one of the determinations of final Kerr parameters in 
Table V. 


mass ratio 

(£,m) 

A 0 (xlO 2 ) tq N m/M 

1:1 

(2,2) 

31.3 ±0.1 11.68 ±0.01 

2:1 

(2,2) 

26.80 ±0.06 11.57 ±0.01 

4:1 

(2,2) 

17.23 ±0.09 11.37 ±0.01 


(2,1) 

4.1 ±0.1 11.42 ±0.11 


(3,3) 

8.2 ±0.2 10.96 ±0.06 


(4,4) 

2.85 ±0.05 10.88 ±0.06 

6:1 

(2,2) 

12.77 ±0.08 11.35 ±0.02 


emitted during the inspiral and merger, there is greater 
uncertainty in the radiation-derived final spin than the 
final mass. For this reason, our preferred measure of 
the final dimensionless spin is the QNM-derived value 
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